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ABSTRACT 


Experiments  indicate  that  a  multigrid-type  cycle  can  be  used  as  an  efficient  precon¬ 
ditioner  in  the  iterative  solution  of  the  discrete  problem  corresponding  to  a  singularly 
perturbed  elliptic  boundary  value  problem.  Motivated  by  a  report  of  Goldstein,  we  ex¬ 
plore  the  theoretical  basis  for  the  efficiency  of  such  a  preconditioner  when  applied  to  a 
model  problem.  The  techniques  developed  are  also  used  to  analyze  a  multigrid  V-cycle 
when  used  alone  as  a  fast  iterative  solver. 


1.  Introduction 


This  work  is  motivated  by  a  report  of  Charles  Goldstein  [7]  in  which  the  author 
discusses  the  task  of  numerically  solving  the  following  elliptic  boundary  value  problem: 


[  u(x)  =  g(x)  on  dft 


+  a0(x)u(x)  =  f(x)  in  ft  C  Ht2 


dii 


(1.1) 


where  x  =  (xi,x2)  €  ft,  0  <  £  <<  1,  the  coefficients  and  data  are  sufficiently  smooth, 
and  cij(x)  >  cq  >  0  in  ft  ,  i  =  0,1,2. 

The  discrete  problem  arising  from  a  typical  discretization  of  (1.1)  on  a  uniform  grid  of 
mesh  size  h,  h  <  e,  is  a  large  system  of  linear  equations.  For  the  solution  of  this  system 
to  approximate  the  solution  of  the  boundary  value  problem  (1.1)  with  a  fixed  accuracy,  we 
must  choose  the  mesh  size  small  for  small  £ ,  specifically,  it  is  sufficient  to  keep  the  ratio 
h/e  fixed  [l],  [11].  In  doing  so,  we  not  only  get  a  much  larger  system,  but  the  resulting 
system  is  also  more  poorly  conditioned. 

With  the  goal  of  trying  to  solve  this  type  of  system,  we  use  the  conjugate  gradient 
algorithm  as  our  iterative  solver.  It  is  known  (e.g.,  [2], [9])  that  if  we  apply  the  method 
of  conjugate  gradients  to  the  problem  Bv  =  F  where  B  is  symmetric,  positive  definite, 
then  the  number  of  iterations,  Nb  ,  required  to  solve  the  system  to  within  a  given  relative 
error,  ||u  —  v'||/llv  “  u°ll  <  V,  is  given  by 

Nb(tj)  <  Cln(2/n)  s/kJB)  (1.2) 

where  I\(B)  =  Amax(B)/Am;n(B) ,  u°  is  the  initial  guess  and  v'  is  the  i  th  approximant 
to  the  solution,  v .  Our  goal  is  to  precondition  the  system  so  that  the  condition  number, 
K  (B*),  of  the  new  system,  BV  =  F’ ,  is  much  smaller  than  K{B)  and  behaves  nicely 
( bounded  or  slowly  increasing)  as  £  and  h  decrease  to  zero. 

It  has  been  observed  experimentally  that  a  certain  multigrid- type  cycle  is  an  inex¬ 
pensive  preconditioner  for  this  system.  The  effectiveness  of  this  preconditioner  is  quite 
sensitive  to  the  choice  of  the  number  of  grids,  fc,  used  in  the  multigrid  process.  Fourier 


analysis  was  used  in  [7]  in  an  attempt  to  prove  that  a  careful  choice  of  the  number  of 
grids  does  guarantee  a  good  preconditioner  in  the  case  where  ft  is  a  rectangle.  Although 
Fourier  analysis  is  routinely  used  to  study  2-grid  multigrid  cycles,  the  k- grid  analysis,  for 
k  >  2  ,  is  quite  unwieldly  and  is  not  usually  attempted.  The  difficulty  arises  from  the 
use  of  coarser  grids  on  which  certain  modes  “alias”  (see  [3])  or  are  “not  visible”  (see  [12]). 
Unfortunately,  this  “aliasing”  was  ignored  in  [7].  The  experimental  evidence  is  so  striking, 
however,  that  it  seemed  worth  trying  to  complete  the  analysis. 

We  examine  the  effectiveness  of  the  multigrid  preconditioner  by  considering  a  special 
case  of  the  boundary  value  problem  (1.1)  with  aj(x)  =  1,  t  sc  0,1,2  ,  &i(x)  =  0,  i  =  1,2, 
ft  =  (0, 1)  x  (0, 1)  and  e  real  and  small.  It  is  for  this  model  operator,  AeL  =  —  s2  A-f/,  that 
we  prove  our  basic  results.  More  general  singularly  perturbed  problems  such  as  variable 
coefficient  and/or  non-symmetric  with  positive  definite  symmetric  part  can  be  analyzed 
using  the  properties  of  the  multigrid  preconditioner  acting  on  A\  together  with  such  ideas 
as  spectral  or  norm  equivalence,  see  [5]  and  [7]. 

Let  h  —  2~n  for  a  positive  integer,  n.  Discretizing  this  model  problem  on  a  uniform 
grid,  ft/,  =  {(/h,  mh )  :  l,m  =  1, 2, . . . ,  2”  —  1}  ,  with  mesh  size,  h,  using  a  standard  5-point 
discretization  of  the  Laplacian  (see  Section  2.1),  we  obtain  the  linear  system 

A\uh  :=  (— e2  Aft  +  I)uh  =  /ft.  (1.3) 

In  Section  3.1  we  define  a  symmetric  linear  operator,  Mjt ,  based  on  multigrid  ideas,  using 
k  —  1  auxiliary  grids  of  larger  mesh  sizes,  2 ph ,  for  p  =  1, 2, . . . ,  k  —  1 .  In  fact,  the  vector 
MkWh  is  essentially  one  “partial”  multigrid  V-cycle  applied  as  if  to  solve  the  problem: 

Ahvh  =  wh,  (1.4) 

starting  with  initial  guess  =  0,  where  A/,  is  the  matrix  resulting  from  the  corresponding 
discretization  of  the  Dirichlet  boundary  value  problem  for  Poisson’s  equation.  In  order  to 
obtain  a  symmetric  operator,  we  take  symmetric  smooths.  I.e.,  if  rp  smooths  are  done 
on  the  p  th  grid  in  the  fine  to  coarse  part  of  the  cycle,  then  rp  smooths  must  be  done  on 
the  p  th  grid  in  the  coarse  to  fine  part.  We  take  a  fixed  rp  =  r  for  all  p  =  0, . . . ,  k  —  1 . 
The  adjective  “partial”  refers  to  the  following  property  of  this  particular  V-cycle:  instead 
of  solving  for  the  coarse  grid  correction  exactly  on  the  coarsest  grid.  2r  iterations  of  the 


2 


smoother  are  applied.  We  choose  the  smoother  to  be  a  damped  Jacobi  iteration  with 
damping  parameter,  cj,  where  0  <  u  <  1.  Taking  u  =  1  would  correspond  to  an 
undamped  Jacobi  iteration,  but  we  exclude  this  choice.  The  choice  u  =  .5  corresponds  to 
a  Richardson  iteration.  Using  Mk  as  a  preconditioner  for  (1.3),  we  claim: 

If  the  mesh  size  on  the  coarsest  grid  is  choosen  to  be  approximately  equal 
to  the  singular  perturbation  parameter,  e,  then  the  condition  number  of 
the  preconditioned  system  is  bounded  independent  of  e  and  h. 

Defining  MjJ  =  M* ,  where  k  is  chosen  so  that  h\  «  e ,  we  justify  this  claim  in  3  steps: 

1.  In  Section  3.2  we  reduce  the  problem  to  finding  appropriate  upper  and  lower 
bounds  for  the  eigenvalues  of  M^Aeh.  Let  q  :  D*  — »  {1,2,  ...,(2n  —  l)2}  : 
(iih,iih)  •-+  qi,i  =  (i’i ,  *2 ) ,  be  a  given  ordering  of  the  (2n  —  l)2  points  of 
fl/, ,  and  let  {oj}  be  a  (given)  complete  set  of  eigenvectors  of  A/,.  Define  a 
(2n  —  l)2  x  (2n  —  l)2  matrix,  M,  by 

=  l1'] 

where 

for  each  i  =  (ti,t2),  j  ~  U1J2)  where  1  <  ii,i2,ji,;2  <  2n  and  (•,•)  is  the 
discrete  -  L2  inner  product.  Using  this  eigenfunction  analysis  (Fourier  analysis), 
the  problem  reduces  to  finding  bounds  on  the  eigenvalues  of  M. .  The  off-diagonal 
elements  of  M  represent  the  “aliasing”. 

2.  In  Section  3.3  we  obtain  a  formula  for  a  bound,  C \  k  ruJ,  such  that,  for  every  i, 

yi  l/^jl  S  Ch,k,r, u/IMmI- 

Therefore  we  have  diagonal  dominance  of  the  matrix,  M. ,  provided  C^,jkir>w, 
where 


can  be  shown  to  be  less  than  one.  The  constant  Ch,k,r,u  is  calculated  for 
r  =  1,2, 3, 4,  u3  =  .5,  .6,  .7,  .8,  .9 ,  h  =  1/2, 1/4, 1/8,. .  .,1/8192  and  all  possi¬ 
ble  corresponding  values  of  k.  All  computed  values  of  Ch,k,r,u  are  less  than  one 
with  the  exception  of  the  case  where  only  one  smoothing  is  used  and  u>  <  .7. 

3.  In  Section  3.5  we  restate  the  bounds  given  in  [7]  on  the  diagonal  entries  of  the 
matrix.  These  bounds  are  used,  combined  with  the  diagonal  dominance,  to  show 
that: 

ci£2  <  Amin(A/£A£)  <  <  c2e2, 

for  constants  Ci ,  c2  >  0 .  The  diagonal  dominance  of  M  is  needed  only  to  guar¬ 
antee  the  positivity  of  the  lower  bound. 

In  Section  4  we  describe  some  experiments  which  illustrate  the  efficiency  of  using  the 
optimal  number  of  grids  in  the  multigrid  preconditioner.  Experimental  comparisons  are 
made  between  three  different  solvers  for  the  model  problem.  In  a  preconditioned  conjugate 
gradient  routine,  two  preconditioners  are  used,  first  the  preconditioner  analyzed  in  this 
paper,  namely  the  preconditioner  based  on  the  Laplacian  with  smoothing  on  the  coarsest 
grid,  and  secondly  a  preconditioner  which  is  based  on  the  model  operator  itself,  solving  on 
the  coarse  grid.  The  third  solver  used  in  the  comparison  is  a  symmetric  multigrid  V-cycie. 

The  techniques  used  in  the  analysis  of  “multigrid-as-a-preconditioner”  can  also  be 
used  to  analyse  “multigrid-as-a-solver” .  This  analysis  is  simpler  than  the  preconditioner 
analysis  since  we  don’t  need  diagonal  dominance  (and  we  don’t  have  it),  see  Section  5.  In 
Section  6  we  show  how  the  k  -grid  convergence  bounds  obtained  in  this  way  compare  to  the 
experimentally  observed  convergence  rates  and  to  V-cycle  convergence  bounds  obtained  by 
other  methods. 


M  WWW  JW«OT«WflW 


2.1 


Notation 

Consider  the  two-dimensional  Dirichlet  problem 

f  -  A  u  =  f  in  ft  =  (0,1)  x  (0,1) 


(2.1) 


1  u  =  0  on  dCl 

where  A  =  d2 /dx2.  We  discretize  this  problem  on  a  family  of  grids.  Let  h  =  2-n, 

as  in  Section  1.  Choose  a  positive  integer  k,  k  <  n.  Define  a  coarse  grid  mesh  size 
/ii  =  2 k~1h  .  In  n  we  define  k  intermediate  grids,  ftp  ,  p  =  1,2 ,...,£  with  mesh  sizes 
hp  =  2 1~phi  .  Clearly  h  =  hk  and 


=  {(*/,  ym )  =  ( lhp ,  mhp )  :  /,  m  -  1, 2, ...,  iVp  -  1}  (2.2) 

where  iVp  =  l//ip  and  p  =  1,2,...,  k  . 

We  define  the  discrete  operator,  Ap ,  which  is  the  negative  of  the  discrete  five  point 
Laplacian,  on  the  grid  using  the  standard  five- point  discretization  of  the  differential 
operator,  —A  (see  e.g.,  [6]).  Each  Ap  is  a  sparse  (Np  —  l)2  x  (iVp  —  l)2  matrix  with  a 
complete  set  of  eigenvectors,  given  by: 

a-p)(m,  n)  =  2sin  (ti7rm/ip)  sin(i27rnfip)  m,  n  =  1, ...,  Np  —  1  .  (2.3) 

where  i  =  (i\,  12)  ,  and  ix,  i2  —  1,2,  ...,iVp  —  1  .  The  corresponding  eigenvalues  are: 

(?)  _  4-2  cos  (iivhp)  -  2  cos  ( i2Trhp ) 

»  L2  #  V 


As  usual,  the  multigrid  operators  we  consider  are  constructed  from  smoothers,  Gp , 
p  =  1, 2, ...,  fc  and  intergrid  transfer  operators,  Ip_ j  and  /p_1  ,  p  =  2, 3, ...,  fc  . 

To  simplify  the  analysis  we  choose  Gp(-,  •)  to  be  a  damped  Jacobi  smoother,  defined 
by 


Gp(uP,  /p)  —  (1  2u>cpAp)up  -f-  2 ujcpfp 

=  Gpup  +  (I  -  Gp)A~l  fp  (2.5) 

where  cp  =  /ip/8,  p  =  1,...,  fc,  and  Gp  is  the  linear  part  of  Gp.  We  require  that 
0  <  u>  <  1.  We  do  not  allow  u;  =  1 ,  which  would  correspond  to  a  Jacobi  iteration.  The 
constant,  cp,  is  approximately  equal  to  the  inverse  of  the  spectral  radius,  p(Ap).  In  fact. 
cpp(Ap)  =  1  -  0 (hp),  and  therefore  Gp  is  a  contraction,  i.e., 


p(I  -  2u'cpAp)  <  1  . 


(2.6) 


We  define  inner  products  and  norms  by: 


and 


for  up,  vp  defined  on  flp  . 


ztttp 


upf  =  (up,up)p, 


(2.7  a) 


(2.76) 


For  the  projection  and  weighting  operators  we  take  Ipx  to  be  linear  interpolation: 


Ip  =  - 
Vi  -  4 


1  2  1 
2  4  2 
J  1  2  1  Up_l 


(2.8a) 


and  Ip  1  to  be  the  adjoint  of  relative  to  the  discrete  -  L 2  inner  products  defined  by 
(2.7a): 

"l  2  1‘ 

2  4  2 
1  2  1 


/p-1  _  J- 
p  16 


P-  1 


(2.86) 


where  we  have  used  the  “distribution”  and  “collection”  stencils  as  in  [10]. 


In  the  eigenfunction  analysis  we  need  some  notation  and  simple  formulas.  Let  i  = 
(11,12).  Define 


and 


dp>  =  cos2 


Jr)  _  ros2  (i2*hp 

**  -cos  ^  — 


(2.9a) 


(2.96) 


A  simple  trigonometric  identity  gives  us 


(p)\ 2 


(2.10a) 


and 


(2.106) 
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The  eigenvalues  of  Ap  can  be  written  as 


v(p)  = 


4(2  -  {<”>  -  ni'*) 


(2.U) 


A  simple  calculation  shows  us  that  the  effect  of  the  projection  on  the  eigenvectors  of 
Ap  can  be  expressed  as(2) 

/'-‘o!"  =  •  (2-12) 

The  corresponding  formulas  for  interpolation  is 


-  +  (1  -  e’xi  - 


(2.13) 


Note  that  eigenvectors  of  Ap  are  also  eigenvectors  of  Gp.  The  eigenvalue,  g[p) ,  of 
Gp,  corresponding  to  ,  is  given  by 


g\p]  =  1  -  2wCpi/,(p), 


where  the  constants  cp  are  related  by 


c0— i  —  4cn. 


(2.14) 


(2.15) 


When  we  apply  the  multigrid  algorithm,  we  transfer  vectors  to  coarser  grids.  In  the 
process  we  lose  information.  In  this  two-dimensional  problem  with  an  (h-2h)  grid  structure 
the  four  (if  ix  ±  Np/2  and  i2  #  Np/ 2)  eigenvectors  a[p) -a\ p^_ti  Ij} ,  -aJp,  ^_tj) 
and  a\pJ,  .  w  _ .  . ,  defined  on  flp ,  are  indistinguishable  on  flp_1  .  There  are  also  2 Np  —  3 
eigenvectors  as  defined  on  f2p  which  are  indistinguishable  from  the  null  vector  as  defined 
on  ftp_1 .  This  phenomenon  is  what  is  referred  to  as  aliasing. 

This  aliasing  plays  an  important  role  in  the  analysis  of  the  multigrid  process  and  we 
introduce  the  following  notation.  Given  two  multi-indices  i  =  ( *"i ,  i"i )  and  j  = 
consider  a\k^  and  a^kK  If  c*.p^  =  ±a(p^  then  we  write  i  ~  j  ( p ).  If  a[p)  and  ai]p)  are  not 
linearly  dependent  then  i  ^  j  ( p ). 

^  In  the  cases  where  |  i  |:=  max(  i\ ,  i2 )  >  \/Np,  one  should  replace  a[p  n  by  its 
proper  (unique)  representation,  aip  ,  where  |  i  ]<  Np-\.  However.  Formula  (2.12)  is 
also  correct  in  this  form. 


I.I  U  LI  U  UH.I 
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2.2  Intergrid  Operator  Identities 

A  multigrid  cycle  consists  of  smoothings  and  intergrid  transfers.  The  smoother  is 
applied  to  reduce  the  high  frequency  (rough)  components  of  the  error.  The  residual  is 
transfered  to  a  coarser  grid  where  solving  exactly  for  the  error  correction  is  less  expensive. 
By  solving  and  then  interpolating  this  coarse  grid  correction  back  to  the  fine  grid,  the  low 
frequency  (smooth)  components  of  the  error  are  reduced.  In  the  boundary  value  problem 
(2.1),  the  eigenfunctions  are  easily  identifiable  as  rough  or  smooth,  being  products  of  sine 
functions.  The  same  is  true  for  the  discrete  operators,  Ap ,  1  <  p  <  k  .  To  gain  insight 
into  the  properties  of  the  multigrid  process  we  study  the  effect  of  a  multigrid  cycle  on  the 
eigenvectors  of  Ak  ■ 

Using  formulas  (2.12)  and  (2.13)  it  is  clear  that  transferring  oqp)  from  flp  to  f lp~l 
and  then  interpolating  back,  results  in  a  linear  combination  of  the  four  eigenvectors  which 
alias  from  Qp  to  f2p_1  .  A  ‘smooth’  eigenvector,  i.e.  and  ^p)  close  to  zero,  picks  up 
‘rougher’  components.  In  the  full  k  -grid  problem  where  there  are  4fc-1  vectors  aliasing 
from  Qk  to  fi1  ,  keeping  track  of  the  aliasing  is  difficult.  Fortunately,  there  are  a  few 
simplifying  features.  The  second  of  the  following  three  Lemmas,  in  particular,  simplifies 
the  analysis.  Define 


v,  =  c+.ch  ■  -  nr' .  i  <  pi  <  P2  <  * 


Lemma  2.1 

If  j  1  (n)  and  )  /  i  (n  +  1)  for  some  0  <  n  <  k  ,  then 


(2.16) 


<«!'> .  ii°T)  = 


Proof  of  Lemma  2.1 


0  if  n  <  p  <  fc; 

m=p+l 


p  <  n 


(2.17) 


Let  j  ~  i  ( n )  and  j  i  ( n  4-  1 )  for  n. ,  0  <  n  <  k. 


For  p  >  n,  the  orthogonality  of  the  a,  gives 


.  )p  =  0. 


For  p  <  n  and  i  /  (0,  0)  ( p ) , 


m=p+l  ^ 

m=p+l  ' 


Since  /;  =  /;/;,  then 


Using  j  ~  i  (n)  and  (2.19)  gives 

'  m=p+l  ' 


If  i  ~  (0,0)  (p)  ,  then  (2.21)  is  trivially  true.  I 


Lemma  2.2 

For  any  n  ,  1  <  n  <  k  ,  and  i  (0,  0)  (n)  , 

£  |(a'n|,/JQf)„|=l. 

(") 

Proof  of  Lemma  2.2 

If  n  =  fc  then  j  ~  i  (1)  implies  j  =  i.  Since  (ajfc\a^)*  =  1,  (2.22)  holds  for 
Assume  £_,~,  (s+i)  I  (<*<,+  1), /*+1a^})J+i  |=  1  for  s,  where  s  <  k. 

i1  =  i  =  (11,12)  , 

?  =  (Na+ 1  —  f  l ,  Z2  )  > 


Define 


Figure  2.1:  A  splitting  of  the  j ,  j  ~  i  (3),  where  s  =  k  —  2. 


The  set  { j  |  j  ~  i  (5)}  can  be  split  into  four  disjoint  subsets  corresponding  to  all 
j  ~  i1  (s  +  1),  j  ~  i2  (3  +  1),  j  ~  i3  (s  4-  1)  and  j  z4  (s  4-  1).  Figure  2.1  shows  this 
schematically  for  the  case  3  =  k  —  2.  Therefore  the  summation  can  be  split  as: 


E  1 1-!*’.  «»!*»>. 


E  1  (2.2-D 

(») 

(  E  +  E  +  E  +  E  ) 

>~i'  (*+l)  J~«3  (s-fl)  j~t3  (j-M)  ;~i4  (s+1) 

ir«!s).4i+1^)4+1 1  • 
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Using  (2.13)  and  the  orthogonality  of  the  of-  ,  the  summation  can  be  written  as: 


E  =  Y  I  (ao*'),I‘k+1a^,),+i  I) 

J~»  (*)  J"*1  (*+l) 

+ (i  -  Y  i  <«!;+,\/;+1<»l‘,>.+.  i 

>~»a  (*+i) 

+  (l-{,“+l,)(l->!.<,+‘))  E  I  (a!.’+1,,/J+IaJ‘)),+1  I 

>~»3  (j+1) 

+ (f.<-+,>)(i  - 1!'+1>)  Y  i  (“l.'+1>.^;+1“!‘,).+i  i  • 

j~t4  (j+1) 

By  the  inductive  hypothesis,  each  summation  on  the  right  hand  side  is  equal  to  one  and 
the  coefficients  also  sum  to  one.  I 


Lemma  2.3 


For  all  n  ,  1  <  n  <  k  ,  and  i  ^  (0, 0)  (n)  , 


(n+l) 


(2.25) 


(n) 

J  +  *  («  +  l) 


Proof  of  Lemma  2.3 

Identity  (2.25)  follows  directly  from  Lemma  2.2  since 

}~ *  («)  J~«  (n) 

(n  +  l) 


-«;”+,V!"+11  E  i<°!n+,’.'rX‘)>. 


;~i  (n+l) 


=  i  -  e 


(n+l)_(n+l) 


«F^-' •* ' '--I'-dv.  .  .  '. 


e 


3.1  Definition  of  the  Preconditioner 


The  multigrid  preconditioner  is  based  on  the  discrete  five  point  Laplacian.  A/*  is  one 
standard  multigrid  symmetric  V-cycle  starting  with  zero  as  the  initial  guess,  except  that 
the  coarse  grid  correction  is  obtained  by  smoothing  instead  of  by  solving  exactly  on  the 
coarsest  grid.  Having  choosen  a  fixed  number  of  grids,  k,  the  multigrid  preconditioner  is 
defined  recursively.  Choose  a  positive  (integer)  number  of  smoothings,  r.  Then  A/*./*  — 
Uk  where  up  (=  Mpfp  ),  for  fp  defined  on  p  =  1, . . . ,  k,  is  given  by: 

1. )  Smooth  r  times  starting  with  initial  guess  =  0: 

Uj>  =  Gp  (0,  fp) .  (3.1a) 

2. )  Compute  the  residual  and  transfer  to  the  coarse  grid: 

rp  =  fp  —  Apup ,  fp—  i  =  Ip~1rp.  (3.1b) 

3. )  Compute  the  coarse  grid  correction: 

If  p  =  2,  up-\  =  ui  =  G]r  (0,  /i )  (3.1c) 

If  p  >  2,  up-i  =  Mp.j/p.j.  (3. Id) 

4. )  Add  the  coarse  grid  correction: 

up  =  up  +  Ipp_xup- 1.  (3.1e) 

5. )  Smooth  r  times  starting  with  initial  guess  =  up: 

up  =  Grp(up,fp).  ( 3.  If ) 

Because  we  have  started  with  an  initial  guess  of  zero,  the  multigrid  preconditioner  is 
a  linear  operator  acting  on  fk  ■  This  definition  of  A/*  can  be  rewritten  as: 

A/p  =  (/  —  G2pr)  Apl  +  GpIp_lMp-iIp~l  Gp  p  =  2,.  ,,k  (3  2) 

A/j  =  (/-Gf^Ar1. 

These  identities  rely  on  the  commutivity  of  Gp  and  Ap.  p  =  1.2,  . k 


and 


3.2  The  Problem 

As  remarked  in  the  introduction,  it  is  sufficient  to  examine  the  effectiveness  of  the 
multigrid  preconditioner  by  considering  the  model  problem  (1.3).  We  take  Q  =  (0.  1)  x 
(0, 1)  and  s  real  and  small.  It  is  for  this  model  operator,  AeL  =  — e2A  +  I ,  that  we  prove 
our  basic  results. 

Define 

Aek  =  £2Ak  +  I.  (3.3) 

Writing  the  symmetric  preconditioner  as  Mk  =  Q*kQk,  the  preconditioned  system  is 
Aeh'v'  =  F'  where  Aeh'  =  QkA\Qk.  Experimental  evidence  suggests  the  following: 

Conjecture: 

Let  r  >  0,  0  <  u>  <  1,  h  >  0  and  e  >  h.  Choose  the  number  of  grid  levels,  k ,  so 
that  hi  =  2 k~1h  m  e.  Define  Mk  =  Mk .  Then  there  exist  constants  ci ,  C2  >  0  such  that 

<  \m*x(MtA\)  <  c252. 

What  we  prove  in  this  paper  is: 

Theorem  3.1 

Let  r  =  1,2, 3, 4  and  u ;  =  .7, .8, .9  or  r  =  2,3,4  and  u  =  .5, .6.  Let  h  >  1/8192  and 
e  >  h.  Choose  k  so  that  h\  =  2 k~1h  ss  e.  Then  there  exist  constants  Ci(h),  c2(h)  >  0 
such  that 


ci(h)e2  <  Amin  (MkAk)  <  Amax  (MkA\)  <  c2(h)e2.  (3.4) 

Remark  3.1 

For  fixed  £ ,  r  and  a/,  numerical  evidence  indicates  that,  as  h  — ►  0. 

ci(h)  — ♦  ci  >  0 

c2(h)  —  c2  >  0. 

Remark  3.2: 

Since  Aeh'  is  similar  to  MkAeh  ,  (3.4)  implies  that  A"(.4^,)  is  bounded  independent 

of  £  . 
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Proof  of  Theorem  3.1: 

Define 

M,;  =  (  Mk  (s2Ak  +  I)  ajfc),  Qji)  )k-  (3.5) 

Because  of  the  aliasing,  Hij  can  be  nonzero  for  j  ^  .  However  if  i  7 C  j  (1)  (i.e.  a(,fc)  and 
atj  '  are  distinguishable  on  the  coarsest  grid)  then  fj.t}  =  0. 

Choose  m  —  where  |mj  :=  max(mi ,  m2)  <  Nk  ■ 

Let  j  1,  J21  •  •  • ,  J4*  —  x  be  some  ordering  of  the  J  '"s''  m(l)  • 

We  now  define  A4m  to  be  a  4fc_1  x  4i_1  matrix  given  by 

(Adm)p?  =  Hjrj,  ■  (3.6) 

We  consider  the  subspaces 

Sm :=  linear  span  a* fc)  :  j  ~  m  (1)})  ,  (3.7) 

where  |m|  <  Nk .  The  Sm  are  orthogonal  (with  respect  to  the  inner  product  defined  by 
(2.7a))  subspaces  and  invariant  under  A/fc(e2A*  +  I).  Therefore  if  we  show  that 

Cl£  5;  ^mia  (^ra)  ^max  (Adm)  ^  C2e  ( 3 . S ) 

for  each  m,  then  (3.4)  will  be  proved. 

By  the  Gershgorin  theorem,  any  eigenvalue.  A,  of  Mm  must  satisfy 

|A-/i„|  <  Y,  M  (3.9) 

;~i  (1) 

;*< 

for  some  i  ~  m  ( 1 ) . 

We  show  that  is  diagonally  row  dominant  and  therefore  we  can  use  information 

about  the  behaviour  of  the  diagonal  entries  of  to  prove  (3  S).  Specifically,  in  Sec¬ 

tion  3  3  we  give  a  computable  formula,  (3.22).  for  a  quantify  C\  k  r  _  .  independent  of  :  . 
such  that 

Y  ^,)\  -  fc.r.»  1  3.10  . 

(1) 


For  certain  choices  of  r  and  u>,  C£  t  p  w  has  been  computed,  for  every  z,  showing  that 
&h,k,r,u>  ■=  sup- C\  k<r  u/  <  1  for  the  k  =  2,  3, . . . ,  12  grid  problems,  using  h  =  2~l  to 
h  =  2-13.  See  Section  3.4.  In  Section  3.5  it  is  shown  that  3  c,  c  >  0  such  that 

cs2  <  min  nn  <  max  an  <  cs2.  (3.11) 

1«I<n*  |>|  <Nk  ~  v  ; 

Combining  (3. 9), (3. 10)  and  (3.11)  we  have,  for  any  eigenvalue,  A,  of  Mm , 

(l  -  Ch<kir,„)  ce2  <  A  <  (l  +  Chtk<rlu,)ce2J  (3.12) 

which  verifies  (3.8)  with  cx  =  (l  -  Ch,k,r,S)  £  and  c2  =  (l  +  Ch,k,r,u)  c. 

Note  that  a  common  factor,  e2v\k)  +  1,  appears  in  all  the  ,  j  ~  z  (1),  therefore 
(3.10)  is  equivalent  to 


(3.13) 


(1) 


Di:= 


(3.14) 


3.3  Bounds  on  the  Off-Diagonal  Elements  of  Mm. 

When  applying  a  multigrid-type  cycle  to  an  eigenvector,  a\k) ,  of  A*,  the  resulting 
vector,  iV/ta|fc),  is  a  linear  combination  of  a[k)  and  all  of  the  other  eigenvectors,  a(k) , 
which  alias  with  on  the  coarsest  grid.  In  this  section  we  give  a  formula  for  a  bound 
on  this  aliasing.  Specifically,  we  find  an  expression,  C\  k  r  w,  where 

J,:=  ]T  |(Mta!k),aJfc))*|  <  a(,fc)  )*.  (3.15) 

(i) 

Let  i  =  (»x,  i2),h,k,r  and  u  be  fixed. 


gp  =  (G;a\p\«\*)p 

(3.16c) 

/2r-l  \ 

(3.16d) 

\  <7=0  / 

and  i/p  =  (ApaSp),aSp))p, 

(3.16e) 

where  the  i,  r,  h  and  u>  dependence  has  been  suppressed  in  the  notation  and  only  the  grid 
level  is  displayed. 

The  following  lemma  gives  a  formula  for  any  entry  in  the  row  of  A4m  corresponding 
to  i,  where  i  ~  m  (1). 


Lemma  3.1 

For  any  j  ~  i  (1), 


(Mka[k\  a{k))k  =  2u jck  ep 

p=i 


(a{p)  Ip  Gr 


p+i  wp+i 


rk~\y~,r  (t)\ 

ik  G*0;  /p- 


(3-17) 


Proof  of  Lemma  3.1 

A  proof  by  induction  shows  that  for  every  s,  2  <  3  <  k , 


(M3a[a\ai;))a=2ucaJ2epl  )  -  (a['\  lpp+lGrp^  ■ 


p=l  \m=p+l 


(3.18) 


Taking  s  =  k  gives  (3.17). 

For  s  =  2,  (3.4),  (3.16)  and  (2.12)  give 

(M2al2) 

+  {{I  -  G\r)  A;' l'2Gr2a[2\  I}Gr2a\2)h 

=  2 ^c2e2(a[2\a*2))2  +  2 — c, e j ^2^2 ^2 ’  /jG^O''21), 


(3  19) 


Substituting  4c2  =  cj  ,  proves  (3.18)  for  s  =  2. 


Assume  (3.18)  is  true  for  s  —  1  grids,  s  >  2.  For  the  s-grid  problem,  (3.4),  (3.16) 
and  (2.12)  give 


<*<'’),  =  {(I-  G]')  A7‘a<'\ 

=  2uc,e,(a\‘>  .a*/’). 


(3.20) 


Using  the  inductive  hypothesis  with  I*~1Graa ^  replacing  and  using  4c,  =  c,_i 

proves  (3.18).  I 


Lemma  3.1  can  be  used  to  get  an  expression  for  Jt,  but  the  summation  over  all 
j  ~  i  (1)  would  be  difficult  to  compute.  Theorem  3.2  shows  that  Jx  can  be  bound  by  an 
expression  which  is  no  more  complicated  than  the  expression  for  D{  —  (JV/fca-^ ,a\k^)k . 

We  claim  that  the  J,  can  be  bounded  by  an  expression  which  is  no  more  complicated 
than  the  expression  for  £>,  (=  ( M kot\k\  a\k^  )k)  : 

Theorem  3.2 


a.)  Dx  =  2 uck  £  ep  ( 

p=  1  V  m=jH-l 


b.  )  J,  <  2 uJCk  £  ep  (  1  -  (  fl  ZrnVm 

p=l  \  \m=p+ 1 


I!  ‘^\9m\£,mrlr 

m  =  p+ 1 


(3.21a) 
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(3.21b 


Proof  of  Theorem  3.2 

a. )  Using  Lemma  3.1  with  j  =  i,  combined  with  equations  (3.16c)  and  (2.12)  proves 

(3.21a). 

b. )  To  prove  (3.21b),  split  the  grid  levels  by  partitioning  the  j  ~  i  (1),  j  i.  See 

Figure  3.1  for  a  schematic  illustration  for  k  =  3.  For  each  n  =  1,2, . . . ,  k  —  1  consider 
the  j ’s  such  that  j  ~  i(n)  but  j  ^  i  (n  -f-  1).  Lemma  3.1,  Lemma  2.1,  Lemma  2.3 
and  (2.6)  lead  to  the  following  bound: 

J,  =  f)  Y.  (3.22) 

n=l  j~i  (n) 


k-l  n  (  n  \  (  k 

—  2u2Cfc  ^  ^  (1  ^n-fl^n+l)  em  n  ^mVm  n 

n=l  p=l  \m=p+l  /  \m=p+l 


Changing  the  order  of  summation  gives 


*— l 

J \  <  2 u;cfc  ^  ep 
P=i 


’k  —  l  n 

^(l-en+lWl)  n  \9m  |£mhm 

.n—p  m=p+ 1 


(3.23) 


Observe  that  the  quantity  in  square  brackets  can  be  simplified  to: 


Figure  3.1:  A  splitting  of  the  j,  j  ~  i  (s),  j  ^  1. 


x  j  ~  1  (1)  (2) 

©  j  ~  i  (2)  ,  j  /  i  (3) 

A  j  ~  i  (3). 


Remark  3.3 

The  constants  Ck,k,r,u  can  now  be  expressed  as 

Ch,k,r,u/  =  SUP  (C/i, jfc,r,w) 


where 


C'h'k.r,  v 


k~  1  /  jt  \  (  k 

E  ep  I  1  —  IT  ^mlm  I  (  II  m  ^?m 

p=l  \  m=p-f  1  J  ymsp+l 

E  Cp  f  n  Wmimnln  ) 

p=l  \m=p+l  / 


Note  that  the  denominator  has  one  more  term  in  the  sum  than  does  the  numerator. 


3.4  Computed  Values  of  the  Off-Diagonal  Bounds 

Ideally,  one  would  like  to  find  analytic  bounds  for  C'h  k  r  w ,  independent  of  i ,  h  and 
k .  On  the  other  hand,  bounds  are  easily  computed  for  any  given  h ,  k ,  r  and  u ;. 

Figures  3. 2-3. 5  indicate  the  dependence  of  on  i  =  (i‘i,i2)  for  A  =  1/64. 

r  =  1 ,  uj  =  .8  and  fc  =2,3,4  and  5  grids.  The  maximum  is  taken  on  the  boundaries  i\  =  1 
or  i2  =  1.  Along  the  boundary  i2  =  1  there  are  2k~ 2  relative  maxima,  for  the  A: -grid 
problem.  (For  all  values  of  h,  k,  r  and  ui  tried,  the  maximum  of  C‘  was  attained  at  (1,  i2) 
and  (i2,  1)  for  some  i2.)  Figures  3. 6-3. 9  show  the  dependence  on  r  for  k  =  4  grids. 

Tables  3. 1-3.8  give  the  calculated  bounds,  sup|j|<1^/l  for  u>  =  .5,  .8  and 

r  =  1 ,2,3  and  4.  The  multi-index  at  which  the  supremum  was  attained  is  listed  below  the 
bound. 

To  find  bounds  for  ui  =  .8  and  r  =  1,2, 3,4,  independent  of  h  and  k ,  we  used 
h  =  1/8192  (which  means  >  67  million  points  on  the  fine  grid).  These  numbers  are 
bounds  for  all  h  >  1/8192  and  all  k  corresponding  to  these  meshsizes.  Observing  the 
asymptotic  behaviour  leads  one  to  believe  that  they  are  also  bounds  for  all  h  <  1/8192 
and  any  number  of  grids,  k.  See  Tables  3.9-3.10. 


Figure  3.2 


Figure  3.3 


Figure  3.4:  C{<k<r<ul  for  h  =  1/64,  r  =  1,  u>  =  .8 


1  SMOOTH 


2  SMOOTHS 


2 


3  SMOOTHS 


4  SMOOTHS 


h 

2  grids 

3  grids 

4  grids 

5  grids 

6  grids 

7  grids 

8  grids 

1/16 

.4640 

(1.9) 

.7165 

(1,11) 

Hi 

1/32 

.4688 

(1,19) 

.7530 

(1,22) 

.9484 

(1,21) 

1.025 

(1,11) 

1/64 

.4707 

(1,37) 

.7632 

(1,45) 

.9942 

(1,41) 

1.149 

(1,21) 

>  1 

1/128 

.4712 

(1,74) 

.7669 

(1,89) 

1.004 

(1,81) 

>  1 

>  1 

>  1 

1/256 

.4712 

(1,149) 

.7669 

(1,178) 

1.006 

(1,163) 

>  1 

>  1 

>  1 

>  1 

1/512 

.4712 

(1,298) 

.7671 

(1,357) 

1.007 

(1,325) 

>  1 

>  1 

>  1 

>  1 

Table  3.2  Ch,k,r,u  cj  =  .5  ,  r  =  2 


h 

2  grids 

3  grids 

4  grids 

5  grids 

6  grids 

7  grids 

8  gTids 

1/16 

1/32 

1/64 

1/128 

1/256 

1/512 

..3115 

(1,9) 

.3196 

(U7) 

.3215 

(1,35) 

.3220 

(1,70) 

,32°1 

(1,139) 

.3221 

(1,278) 

.4296 

(1.5) 

.4574 

(U0) 

.4658 

(U9) 

.4680 

(1,39) 

.4685 

(1,77) 

.4688 

(1,155) 

.4680 

(1,5) 

.5441 

(1,11) 

.5700 

(1,22) 

.5771 

(1,44) 

.5790 
( 1,88) 

.5795 

(1,177) 

.5573 

(1,11) 

.6084 

(1,21) 

.6277 

(1,23) 

.6349 

(1,45) 

.6368 

(1,91) 

.6142 

(1,11) 

.6591 

(1,21) 

.6741 

(1,41) 

.6782 

(1,82) 

.6643 

(1.21) 

.6856 

(1,43) 

.6923 

(1,86) 

.6876 

(1.43) 

.6997 

(1.43) 

Table  3.3  C h,k,r,u>  w  =  . 


2  grids  3  grids  4  grids  5  grids  6  grids  7  grids  8  grids 


1/16 

1/32 

1/64 

1/128 

1/256 

1/512 


.2200 

(1,7) 

.2271 

(1,15) 

.2285 

(1,31) 

.2289 

(1,61) 

.2290 

(1,123) 

.2290 

(1,246) 


.2998 

(1,5) 

.3283 

(1,9) 

.3346 

(1,18) 

.3362 

(1,36) 

.3367 

(1,71) 

.3368 

(1,142) 


(1,5) 

.3694 

(1,10) 

.3756 

(1,20) 

.3773 

(1,39) 

.3777 

(1,78) 


.3961 

(Ml) 

.4093 

(1,22) 

.4130 

(1,44) 

.4139 

(1,88) 


.4165 

(1,21) 

.4243 

(1,23) 

.4279 

(1,46) 


.4170 

(1,21) 

.4297 

(1,21) 

.4356 

(1,41) 


4302 

(1,21) 

.4366 

(1,42) 


2  grids 


1/128 

1/256 

1/512 


niKi 


3  grids 

4  grids 

5  grids 

6  grids 

7  grids 

8  grids 

.2253 

.2368 

(1,3) 

(1,3) 

.2481 

.2659 

.2684 

(1,7) 

(1,5) 

(1,5) 

.2550 

.2823 

.2876 

.2880 

(1,15) 

(1,9) 

(1,10) 

(1,9) 

.2567 

.2868 

.2945 

.3013 

.3016 

(1,31) 

(1,18) 

(1.20) 

(1.11) 

(1,11) 

.2570 

.2880 

2971 

.3081 

3089 

.3090 

(1,63) 

(1,35) 

(1.19) 

(1.23) 

(1.19) 

(  1.19) 

.2571 

.2883 

.2982 

.3099 

3125 

.3131 

(1,126) 

(1,70) 

(1.38) 

(  1.47) 

(1.25) 

( 1.26) 

Table  3.8  Ch,k,r,u>  w  =  .8  ,  r  =  2 


h 

2  grids 

3  grids 

4  grids 

5  grids 

6  grids 

7  grids 

8  grids 

1/16 

.1954 

(1,7) 

.2552 

(1,3) 

.2640 

(1,3) 

1/32 

.2001 

(1,14) 

.2840 

(1,7) 

.2993 

(1,5) 

.3013 

(1,5) 

1/64 

.2013 

(1,28) 

.2932 

(1,15) 

.3223 

(1,9) 

.3252 

(1,10) 

.3268 

(1,9) 

1/128 

.2016 

(1,55) 

.2956 

(1,31) 

.3285 

(1,17) 

.3337 

(1,20) 

.3387 

(Ml) 

.3389 

(1,11) 

1/256 

.2017 

(1,111) 

.2961 

(1,63) 

.3299 

(1,34) 

.3387 

(1,18) 

.3473 

(1,21) 

.3499 

(1,19) 

.3500 

(1,19) 

1/512 

.2018 

(1,222) 

.2963 

(1,127) 

.3304 

(1,69) 

.3402 

(1,36) 

.3497 

(1,42) 

.3537 

(1,39) 

.3538 

(1,38) 

1/8192 

.2018 

.2963 

.3305 

.3407 

.3505 

.3550 

.3581 

Table  3.7  C\  *  r  _  —  =  .8  ,  r  =  3 


A 

2  grids 

3  grids 

4  grids 

-4 

5  grids 

6  grids 

7  grids 

8  grids 

1/16 

■sa 

.1880 

.1891 

ms 

(  1,3) 

(1,3) 

1/32 

.1385 

.1993 

.2077 

.2088 

(1,12) 

(1,7) 

(1,3) 

(1,3) 

1/64 

.1393 

.2032 

.2233 

.2240 

.2243 

(1,24) 

(1,13) 

(1,7) 

(1,7) 

(1.7) 

1/128 

.1395 

.2044 

.2267 

.2296 

.2303 

.2307 

(1,47) 

(1,27) 

(1,15) 

(1.7) 

(1.7) 

(1,7) 

1/256 

.1396 

.2046 

.2278 

.2339 

.2341 

.2353 

.2353 

(1.95) 

(1,54) 

(1,29) 

(1.15) 

(1.14) 

(1,14) 

(1,14) 

1/512 

.1396 

.2047 

.2281 

.2348 

.2357 

.2369 

.2370 

(1,189) 

(1,108) 

(1,58) 

(1.30) 

(U5) 

(1,27) 

(1.29) 

Table  3.8  C\,*,r,w  =  .8  ,  r  =  4 


h 

2  grids 

3  grids 

4  grids 

5  grids 

6  grids 

7  grids 

8  grids 

1/16 

.1482 

1 

ini 

(1,1) 

1 

1/32 

.1057 

.1531 

.1628 

■ 

(1,10) 

(  1.6) 

(1.3) 

1/64 

1067 

.1557 

.1705 

.1716 

.1716 

(1,21) 

(1.12) 

(  1,6) 

(1.6) 

(1.6) 

1/128 

.1068 

.1563 

.1736 

.1761 

.1762 

.1762 

(1.42) 

(  1,24) 

(1.13) 

(  1.7) 

(1.6) 

(1.7) 

1/256 

.1069 

1.56,5 

.1742 

.  1 788 

.1795 

.1797 

.1797 

(1.84) 

(  1.48) 

l  1.26) 

(  1.13) 

(  M3) 

(1.13) 

(1.13) 

1/512 

.1069 

.1565 

.1744 

.1795 

1803 

.1808 

.1810 

(1.168) 

(1.95) 

(1,51) 

(1.27) 

(  1.2.5) 

(1.13) 

I  1.13) 

2  grids 

3  grids 

4  grids 

5  grids 

6  grids 

7  grids 

.3587 

.5196 

.6284 

.6945 

.7487 

.7638 

.2018 

.2963 

.3305 

.3407 

.3505 

.3550 

.1396 

.2047 

.2282 

.2351 

.2370 

.2375 

.1069 

.1566 

.1745 

.1798 

.1812 

.1818 

8  grids 

9  grids 

10  grids 

11  grids 

12  grids 

13  grids 

.7800 

.7896 

.7933 

.7953 

.7948 

.7951 

.3581 

.3589 

.3590 

.3592 

.3592 

.3952 

.2390 

.2394 

.2398 

.2398 

.2398 

.2398 

1821 

.1824 

.1825 

.1825 

.1825 

.1825 

=  12  ,  h  =  1/8192 


ijj  —  6 


i  u  i  ^  ^ u  v  v  u."  v."  h  «.■  I'nn.'W'i.'inj 
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3.5  Bounds  on  the  Diagonal  Elements  of  Mm 

Recall  that  the  diagonal  elements,  p„ ,  of  Mm  where  i  ~  m  (1),  are  given  by, 

H'i  =  (MkAla[k\a[k>)k.  (3.25) 

Since  A\  —  s2  A k  +  I  and  hence 

M»  =  +  l)  A,  (3.26) 

the  bounds  on  the  can  be  obtained  from  suitable  information  about  the  D,  ’s.  The 
following  characterization  of  the  effect  of  the  preconditioner  on  smooth  and  rough  eigen¬ 
vectors  of  Ak  is  central  to  the  analysis  and  was  given  by  Goldstein  in  [7], 


Theorem  3.3 

For  r  >  1 ,  uj  suitably  chosen  and  h  sufficiently  small,  the  Dx ’s  are  positive  real 
numbers  such  that: 


a.) 

A  = 

o  (h\) 

for 

<  d/h\ 

(3.27a) 

b.) 

A  = 

(1  -»l) 
..(fc) 

for 

v\k^  >  d/h\ 

(3.27b) 

where  0  <  p  <  1  and  r?  is  independent  of  h  and  d  is  a  constant. 
We  prove  a  more  explicit  version  of  the  same  result: 

Theorem  3.4 

For  r  >  1 ,  0  <  sjj  <  1  and  a  fixed  constant,  d,  where  j  <  d  <  2, 


a;(  1  ~  uJ)  ^2  ^ 

A 

< 

—h2 

for 

„<*>  <  — 

(3.28a 

a. ) 

max  (2.  d(\  +  ru;))  1  — 

3  1 

*  h] 

b. ) 

>— * 

l 

\ 

A 

< 

1 

for 

„<*>  >  — 

(3.2Sb 

8(1  +  r\jj)v\k) 

*  “  hi  ' 

Proof:  in  appendix. 
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These  theorems  give  us  bounds  on  the  ^lt,  and,  for  example,  Theorem  3.4  leads  to 
the  following  bounds: 


For  <  Sf 


For  i/,(t)  > 


u>(l  —  u)h\ 
max  (2,  d(l  +  ru )) 

du(  1  —  u)e2 
8(1  +  ru) 


2  rud 


S  h..  S 


2  ^1 
£  +7 


2  h\ 

<  Mil  5:  C  +  — 


(3.29a) 


(3.29b) 


Therefore,  taking  h\  =  e,  we  prove  (3.12). 

Using  the  diagonal  dominance  of  the  matrices,  Mm  ,  we  can  estimate  the  dependence 
of  the  condition  number  of  Mk(e2At  +  /)  on  the  ratio  a  =  h\/e 2  from  the  behaviour  of 
the  diagonal  elements,  na.  From  the  inequalities  (3.29)  we  get  an  estimate  for  the  choice 
of  or  which  minimizes  the  condition  number: 


(3.30) 


^optimal  -  8  l  ^ 


This  predicts  that  the  optimal  number  of  grids  decreases  as  the  quantity  ru  increases.  One 
can  also  use  (3.29)  to  show  that  it  is  better  to  choose  too  many  grids,  (a  >  aopt ),  rather 
than  too  few,  (or  <  aopt),  (see  Figure  3.10).  These  observations  all  accurately  describe  the 
experimental  results  —  see  the  next  section. 


approximate 

condition 

number 


a  =  h\/e2 


Figure  3.10:  The  condition  number  estimated  from  the  diagonal  terms 


r 


4.  Multigrid  Preconditioner  —  Experimental  Results. 

Our  numerical  computations  were  carried  out  with  three  objectives  in  mind: 

i)  Observe  the  optimality  of  taking  the  meshsize  on  the  coarsest  grid,  hy  ,  to  approx¬ 
imate  the  singular  perturbation  parameter,  e. 

ii)  Check  the  boundedness  of  the  condition  number  of  the  multigrid-preconditioned 
system  as  e  and  the  fine  grid  meshsize,  h ,  decrease. 

iii)  Compare  the  efficiency  to  other  fast  solvers,  in  particular,  the  corresponding  multi¬ 
grid  algorithm  used  as  an  iterative  solver. 

We  discretize  the  boundary  value  problem: 

(  A‘lu  :=  (—e2  A  +/)u  =  /  in  =  (0, 1)  x  (0, 1) 

1  u  =  0  on  3f2, 

on  a  grid  of  uniform  meshsize,  h ,  as  in  Section  2.1.  Using  the  multigrid  preconditioner,  M k  , 
as  defined  in  Section  3.1,  we  iteratively  solve  the  discrete  problem  using  a  preconditioned 
conjugate  gradient  algorithm.  Recall  that  k  is  the  number  of  grids  used  in  the  multigrid 
algorithm,  hk  =  h,  and  the  smoothers,  Gp  ,  1  <  p  <  fc,used  to  define  Mk ,  depend  on  the 
damping  parameter,  u> ,  and  a  fixed  number  of  smooths  per  iteration,  r .  We  solve 

(e2Ak  +  I)uk  =  (4.2) 

starting  with  initial  guess,  u°k.  We  call  this  iterative  solver  PCCG(-A,sm).  The  “A” 
reminds  us  that  the  multigrid  preconditioner  is  based  on  Ak ,  the  negative  of  the  discrete 
Laplacian,  and  not  on  the  operator  Ak  =  e2  Ak  +  I  and  “sm”  indicates  that  we  smooth 
instead  of  solving  exactly  on  the  coarsest  grid.  Experimentally,  we  find  that  a  reasonably 
good  choice  of  r  and  w  is  r  =  2  and  uj  =  .8  (w  =  .8  is  optimal  for  the  corresponding 
2-grid  multigrid  solver,  see  [12]). 

We  first  consider  solving  (4.2)  with  Fk  =  1.  For  h  =  1/64  we  show  the  dependence 
of  the  number  of  iterations  required  to  reduce  the  norm  of  the  residual  by  a  factor  of  10-6 
on  the  choice  of  e  and  h j .  See  Table  4.1  .  For  given  e  and  h,  the  number  of  iterations 
listed  is  the  largest  observed  for  various  choices  of  u°k  .  Note,  in  particular,  the  cases  where 
hi  =  e. 

Table  4.2  displays  the  number  of  iterations  required  to  reduce  the  relative  error  by  a 
factor  of  10-6  for  various  choices  of  h  and  e,  taking  hA  =  s.  Here  we  used  Fk  = 


.  *  *■*  .  *.  i  -  *■  _  •  *-  _  *  _  1.  _  p* 
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Finally,  we  compare  the  efficiency  of  PCCG(— A,sm)  to  other  elliptic  solvers.  We 
take  h  =  1/64,  e  =  1/8,  Fk  =  1  and  an  initial  guess  consisting  of  a  smooth  and  a  rough 
component,  namely: 

u°k  =  10  +  20  cos(647rx)  cos(647ry). 

We  consider  a  symmetric  V-cycle,  which  is  a  fast  iterative  solver  for  equation  (4.1),  where 
we  solve  exactly  on  the  coarsest  grid  (we  use  a  symmetric  band  solver  to  invert  £2  A\  +  I  ). 
We  denote  this  algorithm  by  MULT.  For  comparison,  an  (extreme)  choice  of  a  precondi¬ 
tioner  for  the  preconditioned  conjugate  gradient  algorithm  is  considered,  where  the  pre- 
conditioner  is  based  on  A\  instead  of  Ak  and  we  solve  exactly  on  the  coarsest  grid.  In 
other  words,  this  preconditioner  consists  of  one  cycle  of  the  solver,  MULT,  starting  with 
initial  guess  of  zero.  This  algorithm  is  called  PCCG(-e2  A  +  I , so).  Of  course  we  expect 
the  behaviour  of  this  preconditioner  to  be  better  than  that  of  the  simpler  ( —  A,sm)  pre- 
conditioner,  but  we  have  the  added  expense  of  a  coarse  grid  solve  and  (slightly)  more 
complicated  operator.  Of  interest  to  us  here  is  that  PCCG(  — e2  A  -fi  I, so)  is  not  a  signif¬ 
icant  improvement  over  PCCG(  —  A,sm)  if  the  optimal  choice  of  the  number  of  grids  is 
used. 

In  a  conjugate  gradient  algorithm,  the  error  reduction  factor,  ||ejfc||/||e*-i  |] ,  typically 
decreases  as  k  increases,  whereas  for  a  multigrid  algorithm  the  error  reduction  factor 
increases  as  k  increases.  Therefore  the  preconditioned  conjugate  gradient  routines  will 
be  more  competitive  when  a  large  reduction  in  the  relative  residual  is  required  and  the 
multigrid  algorithm  is  more  competitive  when  a  smaller  reduction  in  the  relative  residual 
is  required. 

We  also  observe  that  increasing  the  number  of  smoothings  per  grid  level  will  im¬ 
prove  the  performance  of  MULT  more  than  it  will  improve  the  performance  of  the 
PCCG(  —  A,sm)  algorithm.  Similarly,  optimizing  the  choice  of  the  damping  parameter, 
w,  will  improve  MULT  more  than  it  will  improve  PCCG(  —  A,sm). 

Furthermore,  one  should  keep  in  mind  that,  though  it  is  difficult  to  improve  the 
behaviour  of  the  multigrid  preconditioner,  it  is  quite  obvious  how  to  improve  the  multigrid 
solver.  Using  better  smoothers,  or  using  a  full  multigrid  algorithm  (FMG)  will  dramatically 
improve  the  convergence  rate. 

Our  first  comparison  is  made  with  parameters  which  should  give  the  PCCG(  -A.sm) 


algorithm  an  advantage.  We  therefore  consider  a  relatively  inefficient  choice  of  the  damping 
parameter,  u>  —  .5,  and  require  the  norm  of  the  residual  to  be  reduced  by  a  factor  of  10-12  . 
The  total  cpu  time  (seconds)  is  recorded  in  Figure  4.3,  with  the  number  of  iterations  given 
in  parentheses  next  to  the  time.  The  PCCG(-A,sm)  algorithm  appears  to  be  competitive 
with  MULT,  at  least  for  this  meshsize,  h.  The  PCCG(-e2A  ■+  /,so)  algorithm  is  only 
slightly  faster. 

We  then  take  a  more  reasonable  value  of  u>  =  .8  and  require  the  norm  of  the  residua) 
to  be  reduced  by  a  factor  of  10-6 .  The  total  cpu  time  is  recorded  in  Figure  4.4.  The 
multigrid  solver,  MULT,  is  now  the  best  choice. 

All  computations  were  done  on  a  VAX  11/780. 

We  end  this  section  with  a  few  comments  on  the  choice  of  using  multigrid  by  itself  as 
a  solver,  or  using  multigrid  (based  on  a  simpler  operator)  as  a  preconditioner: 

-  For  the  model  problem  (8.1),  our  experiments  indicate  that,  for  modest  values  of  h 
and  e ,  a  good  multigrid  algorithm  is  more  efficient  than  a  multigrid-preconditioned 
conjugate  gradient  algorithm. 

-  In  a  true  variable  coefficient  problem,  (1.1),  the  multigrid  preconditioner  has  the  ad¬ 
vantage  of  being  based  on  a  constant  coefficent  operator.  In  this  case,  using  multigrid 
as  a  preconditioner  should  be  more  competitive  than  in  the  model  problem  case.  It 
is  doubtful  whether  the  multigrid  preconditioner  could  outperform  a  good  multigrid 
solver  even  in  this  case,  but  more  testing  would  need  to  be  done. 

-  In  an  indefinite  problem,  where  multigrid  solvers  are  more  troublesome,  one  of  the 
preconditioned  conjugate  gradient  routines  for  indefinite  problems  might  be  preferable. 


Table  4.1  Optimality  of  choosing  h\  sa  e. 

Largest  (observed)  #  of  iterations  required  for  llrfcl!/llrfcll  <  io-6. 

Fjt  =  1  ,  u  =  .8  ,  r  —  2 


e  =  1/2 

£  =  1/4 

£  =  1/8 

pisi 

>  20 

>  20 

20 

ESS 

12 

12 

10 

9 

8 

8 

BS 

7 

7 

9 

mm 

7 

8 

9 

Table  4.2  Boundedness  of  condition  number  independent  of  h  and  e  taking 
Largest  (observed)  #  of  iterations  required  for  \\uk  —  u*fe||/||ujfc  —  u°k\\  <  10-6. 

Fk  =  0  ,  u  =  .8  ,  r  -  2 


h 

£  =  1/4 

00 

r-4 

II 

£  =  1/16 

£  =  1/32 

5 

6 

■ 

6 

6 

6 

■ 

6 

6 

6 

6 

Table  4.3  Experimental  comparisons  of  approximate  cpu  time  (sec). 
Approximate  cpu  time  (no.  of  iterations)  required  for  j|resfc||/j|reso||  <  10~12. 


Fk  =  1  ,  w  =  .5  ,  r  =  2 

e  =  l/8  ,  /i  =  1/64  ,  u°k  =  10  4-  20  ■  cos  64m  cos  647ry 


#  of  grids 

MULT:V(2,2) 

PCCG(- A,sm) 

PCCG(-£2A  +  l,  so) 

2 

61.3  (20) 

-  (>20) 

53. 4  (io) 

4 

44.2  (2x) 

40.6  (n) 

39.2  (10) 

6 

44.4  (M) 

44.8  {12) 

39.5  (10) 

Table  4.4  Experimental  comparisons  of  approximate  cpu  time  (sec). 
Approximate  cpu  time  (no.  of  iterations)  required  for  ||rest||/||res0||  <  10~6. 

Fk  =  1  ,  w  =  .8  r  =  2 

e  =  1/8  ,  h  =  1/64  ,  u°k  =  10  +  20  •  cos  64xx  cos 


#  of  grids 

MULT:V(2,2) 

PCCG(-A.sm) 

PCCG(-£*’A  +  /  .so: 

5.1  V-cycle  Convergence  Bounds 

In  this  section  we  briefly  describe  the  results  of  applying  the  same  techniques,  in 
particular  Lemma  2.2,  to  obtain  bounds  on  the  asymptotic  convergence  rates  for  multigrid 
V-cycles  used  to  solve  the  Dirichlet  problem  for  Poisson’s  equation  in  the  unit  square.  The 
analysis  is  simpler  in  this  case  because  we  don’t  need  diagonal  dominance.  Instead,  we 
numerically  evaluate  the  ||  •  norm  of  the  appropriate  matrix  (i.e.,  the  largest  row  sum 
of  absolute  values)  which  is  a  bound  on  the  spectral  radius.  We  present  the  details  of  this 
analysis  in  Section  5.2.  We  first  define  our  basic  multigrid  V-cycle  applied  to  the  linear 
system 

BkUk  —  Fk  (5.1) 

starting  with  initial  guess,  u°k,  with  auxiliary  problems,  BpUp  =  /p,  p  =  1,2, . . . ,  k  —  1, 
corresponding  to  discretizations  on  the  coarser  gnds. 

1.  Initialize: 

fk  <-  Fk 

o 

Ufc  <-  uk 

2.  Update: 

Uk  uk 

where  each  up ,  p  =  2, 3, . . . ,  k  is  defined  recursively  by: 

(a.)  Smooth  r  times  starting  with  initial  guess  =  up: 

up  —  Grp(up,  fp)  (5.2 a) 

(b.)  Compute  the  residual  and  transfer  to  the  next  coarser  grid: 

Tp  —  fp  —  Bpup,  fp-1  =  Ip~lrp  (5.26) 

(c.)  If  p  >  2  then  return  to  (a.)  to  evaluate  up_j .  If  p  =  2  then: 

ui  =  B~lh  (5.2  c) 


(d.)  Add  the  coarse  grid  correction: 

Up  =  Up  Ip_]Up—  1 


(5.2  d) 
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(e.)  Smooth  s  times  starting  with  initial  guess  =  up: 

Up  —  Gp(up,  fp) 


(5.2c) 


For  the  model  problem  analysis,  we  take  ftp,  Ap,  Ip_x,  1 p  1  and  Gp  as  defined  in 
Section  2.1. 

5.2  Error  Analysis 

Bounds  on  the  asymptotic  convergence  factors  of  the  multigrid  cycles  M\h,k  r,w  can 
be  found  in  the  following  manner.  Let  ejt  =  Uk  —  Uk  be  the  initial  error  and  ik  —  Uk  -  uk 
be  the  error  after  one  multigrid  cycle,  where  £7*  satisfies  AkUk  =  fk  ■  In  terms  of  the 
errors,  definition  (5.2)  becomes: 


( a)  For  p  =  k  ,k  —  1 ,  •  •  •  ,2 


(b)  For  p=l 


(c)  For  p  —  2,-  ■  ■  ,k 


£p  —  Gp£p  i 
o-l  =  Ap£t 


e,  =  0. 


ep  —  ep  7p_j(ep_i  £p- 1) 


Recall  that  Gp  is  the  linear  part  of  Gp.  If  Mk£k  =  £k  ,  then  Mk  is  defined  recursively 


Mp  =  Grp-  i;_x(I  -  Mp~l)A~_1I£~1GTpAp  ,  2  <  p  <  k 

M1  =  0  . 


Note  that  the'a^  are  eigenvectors  of  Ak  and  Gk  ,  but  not  of  Mk  .  Define 


(5.3a) 
(5.3  b) 


S,  =  linearspan  {q^  :  j  ~  i  (1)}  .  (5  4) 

By  formulas  (2  12)  and  (2.13)  we  see  that  the  S,  are  orthogonal  subspaces  which  are 
invariant  under  Mk  .  Therefore  a  basis  of  eigenvectors,  ,  of  AIk  exists  such  that. 


each  Up  can  be  written  as 


•  (i) 
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for  some  i,  |  i  |<  Nk  ,  where  a e  1R.  Since  the  a,  ;  are  orthonormal  with  respect  to  the 
discrete  L2  inner  product,  then 

(.W‘o„»,)4=(  V  a^M'af  ,  £  )» 

j~i  (1)  m~«  (1) 

=  I]  ^  a]lx(Mka(-\a(^)k  =  AM  ^  (5.6) 

m~i  (1)  (1)  n~)  (1) 

where  AM  is  the  eigenvalue  of  Mk  corresponding  to  . 

A  bound  on  the  AM ’s  will  be  a  bound  on  the  asymptotic  convergence  rate  of  the 
multigrid  cycle.  Let  M,  be  the  4fc_1  x  4*-1  matrix  with  (Ad,)Pi?  =  a**')*  with 

some  ordering  of  all  the  j  ~  t  (1). 

Remark  5.1  Note  that  for  some  i’s,  these  jp’s  are  not  necessarily  unique.  For  example, 
if  i  =  (Nk/ 2, 1)  then  (JVfc/2, 1)  =  (JV*  -  Nk/ 2, 1)  . 


Remark  5.2  The  diagonal  elements  of  A are  the  Rayleigh  quotients, 

<«*■«*>* 

and  the  off-diagonal  elements  are  the  contribution  from  the  aliasing  vectors. 
By  Gershgorin’s  theorem,  any  eigenvalue  A  of  M,  must  satisfy 


i  A  -  (il/*aS,*», «<,*>)*  |<  £  I  I  (5.7) 

>~n  (1) 


for  some  n  ~  i  (1)  .  Therefore  a  bound  on  the  asymptotic  convergence  rate,  p  ,  is  given 
by 


p  <  max 
l-KA'i 


=  max 
I » I  <  jv* 


£  i  i 

>~«  0) 


(5.8) 


for  the  fc-grid  problem  with  meshsize  hk  =  l/Nk  on  the  fine  grid. 

In  section  5.3  we  derive  formulas  for  a  bound  on  the  righthand  side  of  (5.8). 


5.3  Derivation  of  Bounds  on  the  Convergence  Rate 

For  a  fixed  fine  meshsize  h ,  a  given  number  of  grids  k ,  r  smoothings  and  a  damped 
Jacobi  parameter  u>,  we  derive  formulas  for  a  constant  Ck,r,h,u  <  1)  independent  of  i 
which  is  a  bound  on  the  asymptotic  convergence  rate.  In  Section  5.4  we  give  values  of 
these  constants  for  various  values  of  h,  k  and  r  using  a  typical  value  of  ui . 

By  (5.8)  it  is  enough  to  bound  ^  |  (Mka\k\a^)k  j  independent  of  i  .  Divide 


the  sum  into  two  parts, 


;~«  (i) 


;~i  (1) 


+  E  i  («*<.!*’ . 


J~«  (0) 

='•  D{  +  J{  , 


where  D ,  is  the  “diagonal  part”  and  Ji  is  the  “aliasing  part”  of  the  sum. 
Let  i  =  ( i  j ,  i2 ) ,  k ,  r ,  h  and  u  be  fixed.  Define 


f,  =  <.<rt  =  cos2 


(5.10a) 


(p)  2  (l 27r^1 

Ip  =  =  COS  — — - 


^  /▼— n  ' 


(5.106) 


(5.10c) 


(5.10  d) 


vp  ~  V\P)  > 


[5.10c) 


where  the  i,  r,  h  and  u>  dependence  has  been  suppressed  m  the  notation  and  only  the 
grid  level  is  displayed. 

We  have  the  following  theorem. 


Theorem  5.1 


/  *  \ 

D,  —  gk~ 2u:ckvk  ^  *  fpi  Jj[  4pm^mr7mj  — 

p= 2  ^  m  =  p-t- 1  / 


CkVk 


P+1 


Clt'l 


II 


m  =  2 


and 


*=p+ 


p+i 


+  ~  1  1  —  (  TT  im^m 

C\V\  V  A  A 


m= 2 


'5  11a 


k  \  /  * 

J,  <  2ucki/k  EM>-(  n  )  (  1 1  4  |  Qm  I  l?m 

p=2  '  m=p-fl  ^  m  =  p-fl 


(5.116) 


(  1 1  4  |  gm  |  ifm^r 
'  m=2 


Remark  5.3  Theorem  5.1  allows  us  to  obtain  a  bound  on  the  asymptotic  convergence 
rate  that  is  no  more  complicated  than  the  diagonal  elements  themselves. 

Before  proving  Theorem  5.1  we  find  expressions  for  the  inner  products 

/  ( k)\ 

{M  a,  ,  a j  }k- 


Lemma  5.1 

For  any  j  ~  i  ( 1 ) , 


*-i 


[XIka\k\a(k))k  =  gk{a(k\a(}k))k  -  2ujckvk  /p  f  fj  4gm^mJ7m  J  (a(,p) . 

p= 2  m  =  p+l  ' 

k 


Ckl'k 


(5.12a) 


Proof  of  Lemma  5.1 


We  prove  by  induction  that  for  every  s  <  k  , 


(M3a[s)  ,a}3))3  =  g3(a[a) ,  a^)  s  —2ujc,vs  £  fJ  4  gmUVm)  (*\P' .  I*a\")p 

p=2  ^  m  =  p+ 1  ' 

~~~  (  II  (<*[ '  ' .  j\  a' * '  )  i  (5.126) 

m  1 


Taking  s  =  k  gives  (5.12a). 

We  start  with  $  —  2  From  (5.3),  (5.10)  and  (2.12), 


,,, 2  (2)  (2), 

(A/  q,  ,  a  ): 


/  sir  (  1 )  ( 2)  \  /j—  1  rMr  j  (2)  rl  ( 2) \ 

=  (G2a,  ,  )2  -  (.4j  2 2  J2a;  )  i 


(2) 


f2h 


(5.13) 


_  „  /  (2)  (2K  ^2  ,  ,  (1)  rl  (2)  \ 

=  9 2\a,  ,a;  )2  -  —  gibni  {a,  ,/2Qr;  )i  • 


Using  4c2  =  C!  gives  us  (5  12a)  for  k  =  2. 

Assume  (5.12a)  is  true  for  k  =  s  —  1  grids,  3  >  3.  For  the  s-grid  problem,  (5.3), 
(5.10)  and  (2.12)  give 


(A Tal’X’*),  =  (G:a\’\a\”),  -  ((/ - 


.(»)  <*) 


1-1  r*-i< 


(») 


>~i  (*)\ 

Q  , 


‘'j  c  /  (*-i)  r»-l 

—  isHtSsiot,  ,  Is  Otj  )s  —  1 

"j-l 


+ 


—^ns9s(Ms' 

V3-1 


a 


(3-1) 

1  ’ 


ir 


(5.14) 


We  factor  1  -  g,-\  =  2^c,_i  i/,_  { /,_]  .  Using  the  inductive  hypothesis  and  using  4cs  = 
c,_i  finishes  the  proof.  I 


Proof  of  Theorem  5.1 

(a)  Using  Lemma  5  1  with  j  =  i.(5  10c)  and  (2  12)  proves  ( 5  11a) 

(b)  To  prove  (5  lib),  split  the  grid  levels  bv  partitioning  the  j  ~  1  (1),  j  /  /.  See 

Figure  3.1  for  a  schematic  illustration  for  k  =  3  For  each  n  =  1 .  .6-1  consider 

the  j  s  such  that  y  ~  /  (ni  but  j  /  /  t  n  1  )  Lemma  51.  Lemma  2  1  and  Lemma 
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2.3  show  that 


*-i 


J.-Z  Z  k^X*’). 


(5.15) 


n=l  j~i (n) 
X*  (n  +  1) 


k  —  1  n  y  n  \  /  * 

<  2u>ckvk  X/1  sn+l^n+l)  X  'M  n  £"»'?"»)(  n  I  9m  |  s m  T]m 

n=  1  p—  1  m  =  p-fl  m=p+ 1 

»  fc  —  1  ✓  n  \1/^ 

^(1  -  tn+\Tln+\)(  JJ  £mbn  J  (  1 1  4  I  !?m  i  r?: 

•  n=l  ^  m=l  m=l 


cm 

Cl  I'l 


Changing  the  order  of  summation  gives 


fc-i  k- i  /  n  \  I  /  * 

<  2ujckuk  X  X^1  ”  l'n+1M  n  n  4  I  5-m  |  6^, 

p=2  n=p  '  m  =  p-f  1  '  ^  ^  m=p-M 

*  /  n  \  1  /  ^ 

5^(1  -  Cn+ir?n+l)f  J| 

■  n=l  ^  m=p+l  '  ^01=1 


(5.16) 


|  ckuk 

Cl  1/1 

The  quantitites  in  the  square  brackets  in  (5.16)  equal 

k 

1  —  ||  (m’Jm 


and 

* 

1  |j[  £mbm 

mal 

respectively,  and  therefore  (5.11b)  has  been  proved.  I 

We  use  this  theorem  to  find  bounds  on  the  multigrid  V-cycle  asymptotic  convergence 
rate  for  the  /t-grid  problem  with  a  given  damped  Jacobi  parameter  u)  and  r  iterations  per 
smooth.  The  results  are  given  in  the  next  section. 


5.4  Computed  values  of  the  asymptotic  convergence  bounds 

Ideally,  one  would  be  able  to  compute  A: -grid  convergence  bounds  independent  of 
h  .  The  4fc_1  x  4*'1  matrix,  At,,  can  be  written  as  a  4*-1  x  4*;-1  matrix.  £.17). 
with  variable  entries  depending  on  the  continuous  variables  £  and  17  f  (0.  1)  evaluated  at 


In  the  two  grid  case  one  could  get  an  analytic  formula  for  the  characteristic  equation 
of  (a  polynomial  of  degree  4  for  fixed  £,  n  j.  find  analytic  expressions  for  the 

eigenvalues  and  then  find  the  supreraum  of  these  expressions  over  all  £  and  17  in  the  unit 
square.  This  would  give  an  exact  2-grid  asymptotic  convergence  rate  independent  of  h  .  In 
practice  this  is  too  much  work  even  in  the  simple  2-grid  case.  Instead,  one  chooses  a  value 
of  h  —  l/N  and  computes  the  spectral  radii  of  .Vf,  for  each  1,  |  t  |<  „V ,  keeping  track 
of  the  largest.  One  then  repeats  the  procedure  for  different  values  of  h  and  so  constructs 
a  table  as  in  [12]  see  Table  5.1  From  such  tables  one  can  predict  the  h -independent 
convergence  rates. 

In  the  k  grid  problem,  k  >  2  each  M,  is  a  4t_1  x  4fc_1  matrix  and  therefore 
computing  the  spectial  radius  for  each  i,  |  i  |  <  IV  is  expensive,  especially  for  small 
h  .  We  therefore  use  Theorem  5.1  and  Gershgorin’s  Theorem  to  compute  a  bound  on 
the  spectral  radius  of  for  each  1  .  This  amounts  to  roughly  twice  the  work  of  just 
evaluating  the  diagonal  elements. 

The  sharpest  bounds  on  the  asymptotic  convergence  rates  for  the  analysis  of  the 
V -cycle  are  obtained  by  these  techniques  when  no  smoothing  is  performed  on  the  coarse- 
to-fine  part  of  the  cycle,  i.e.,  s  =  0  in  step  d.  This  is  called  an  \I\  cycle.  The  symmetric 
cycle,  i.e.,  s  =  r,  is  called  an  MG  cycle.  We  consider  two  discretizations  of  the  Laplacian. 
the  five  point  discretization,  Bp  =  Ap ,  as  given  in  Section  2,  and  a  certain  nine  point 
discretization  given  by  the  following  stencil: 

(5.17) 

The  corresponding  V-cycles  will  be  denoted  by,  e.g.,  A/5\,  or  MG$ ,  to  indicated  which 
discretization  is  being  used. 

We  consider  a  M$\  algorithm  and  compare  our  theoretical  bounds  to  the  experimen¬ 
tally  observed  asymptotic  convergence  rates.  In  order  to  compare  our  two  grid  bounds 
to  the  exact  two  grid  convergence  rates  obtained  by  the  model  problem  analysis  in  [8], 
we  consider  a  damped  Jacobi  paramete  —  4/5.  Experimentally,  this  is  a  good  choice, 
though  its  optimality  depends  on  the  number  of  smoothings  and  the  number  of  grids  We 
take  r  =  1 ,  2.  3  or  4  smoothings  (smoothing  only  from  fine  to  coarse  meshes).  Tables  5.2- 
5.5  show  the  convergence  bounds  for  commonly  used  meshsizes.  Table  5.6  indicates  the 


limiting  behaviour  of  these  rates  for  very  small  h  and  large  number  of  grids  The  experi¬ 
mentally  observed  asymptotic  convergence  rates  are  shown  in  Table  5.7  for  r  =  1.2.  3.  4. 
_■  =  4/5  and  h  =  1/64.  For  exact  two  grid  convergence  rates,  see  Table  5  1 

In  practice,  as  k  increases  there  is  not  as  much  degredation  in  the  convergence  rate 
as  Tables  5. 1-5. 7  would  indicate. 

We  compare  our  bounds  to  the  finite  element  bounds  of  [8],  using  the  A/Gg  cycle 
given  by  taking  Bp  =  Ap  and  s  =  r.  The  comparison  is  possible  because  the  operators 


-4P  satisfy: 


for  p=  1.2. 


(5.18) 


'  F  *  ~  p  Fp—  i  r  ‘ ~ v  -*w' 

Eigenvectors  of  ,4P  are  also  eigenvectors  of  Ap.  We  also  note  that  for  a  symmetric  V- 
cycle,  convergence  bounds  in  the  energy  norm  are  equivalent  to  asymptotic  convergence 
bounds  given  by  the  spectral  radius.  Our  bounds  are  given  in  Table  5.8  for  u  =  3/4. 
h  =  1/64.  and  r  =  1,2.  3,  4.  In  the  next  to  the  last  column  of  Table  5.8  we  show  the 
bounds  (which  are  independent  of  the  number  of  grids  used)  obtained  by  the  methods  of 
[8].  We  also  calculate  the  exact  two  grid  convergence  rates  for  A/Gg,  as  in  [12].  These 
numbers  are  given  in  the  last  column  of  Table  5.8.  In  this  symmetric  case,  at  least  for  small 
r ,  our  bounds  are  larger  than  the  finite  element  bounds  because  in  the  Fourier  analysis  we 
essentially  throw  away  the  post  smoothing  factors  in  the  off-diagonal  terms  in  order  to  be 
able  to  apply  Lemma  5.1. 


Table  5.1  AT, \  Two  grid  asymptotic  convergence  rates  u>  =  .8 


.592 

.351 

.598 

.358 

.600 

.359 

.600 

.360 

)le  5.5  M$\  Asymptotic  convergence  bounds  u  =  .8  ,  r  =  4 


Table  5.7  A/sV  Experimental  asymptotic  convergence  rates 


U>  =  .8,  h  =  1/64 


2  grids 

3  grids 

4  grids 

5  grids 

6  grids 

.600 

.600 

.600 

.600 

.600 

.360 

.360 

.360 

.360 

.360 

.216 

.228 

.233 

.242 

.246 

.137 

.158 

.171 

.181 

.193 

Table  5.8 


A  comparison  of  the  theoretical  bounds 


w  =  .75,  h  =  1/64 


2  grids 

3  grids 

4  grids 

5  grids 

6  grids 

.686 

.717 

.816 

.860 

.879 

.275 

.299 

.348 

.362 

.364 

.121 

.147 

.161 

.162 

.162 

.079 

.114 

.124 

.124 

.124 

bounds 
from  [8] 


exact  2  grid 
conv.  rates 


I 


APPENDIX 


A.l  Proof  of  Theorem  3.4 

Fix  i  =  (ij,  i2),h,k,r  and  u>  as  in  Section  3.3.  Define  £p,tjp,  gp,ep  and  up  as  is 
(3.16a-e).  As  seen  in  the  proof  of  Lemma  3.1, 


£)<'":=  {Mpa^\a^)r  =  {(I-C*;)A-'c/?\aY>), 

•+  (Mr.yqr'G',a'f\ /'-‘CJa  J'’),.,. 


(A.l) 


Therefore  a  recursion  formula  for  D\p}  is 
D\p)  =  ap  +  bpD\p-l) 
where  ap  and  bp  are  given  by: 


P=  1,2, .  ..,k 


(A. 2) 


ap  =  2tjjcpep 

P  =  1 ,...,k 

(A. 3a) 

P  = 

(A. 3b) 

bo  =  0. 

(A. 3c) 

The  following  four  lemmas  are 

all  proved  by  direct  calculation. 

Lemma  A.l 

For  each  p  —  1, 2, . . . ,  k 

a.)  ap  <  Aruicp 

(A. 4a) 

and  b.)  ap  >  4cj(1  —  u>)cp. 

(A. 4b) 

Lemma  A. 2 

For  each  p  =  2, . . . ,  k 


a.)  bp  <  1 


(A. 5a) 


and  if  cpvp  <  1/4, 

b.)  bp  >  (1  -  4(1  +  ru)cpvp) 

Lemma  A. 3 

For  each  p  =  2, . . . ,  k 


a.)  vpjvp. _!  < 


tpVp 


and  b.)  vpjvp-X  >  1. 


Lemma  A. 4 


0  0  1 

If  4^+r  -  cvvp  -  P  ~  2>  •  •  •  1  k  and  -  <  0  <  2,  then 


a.)  Cp_nt/p_n  < 
and  b.)  cp_n j/p_n  > 


0 


4a— n 

$ 


1-?.  ' 


4a-n+l  V  3  40,-71+1  i  • 


Proof  of  Lemma  A.l 

Inequality  (A. 4a)  follows  immediately  from  the  inequality 

1  —  (1  —  x)2r  <  2 rx 

since  |1  —  x\  <  1  where  x  =  2ucpvp. 

Using  the  inequality 

1  —  (1  —  x)2r  >  x{2  —  x)  for  all  x  such  that  1 1  —  x |  <  1, 

it  is  clear  that 

ap  >  2+cp  (2  -  2ucpup) , 

from  which  follows  (A. 4b)  since  0  <  cpvp  <1.  I 


(A. 5b) 

(A. 6a) 
(A. 6b) 

(A. 7a) 
(A. 7b) 

(A. 8) 

(A. 9) 

(A. 10) 
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Proof  of  Lemma  A. 2 


Since 


CpVp  = 


_  "  £p  V  p 


where  0  <  £p,  r?p  <  1  and  |gp|  =  |1  —  2u>cvvp\  <  1,  (A. 5a)  is  obvious. 
If  CpUp  <  1/4,  then  (1  —  £p)  and  (1  —  r?p)  <  1/2  and  therefore 

/„(p-i)  Jp-D\  _  i 

[Otj  ,  Of,  /p— 1  t. 

Moreover, 

tlvl  >  1  ~  4cpi/p. 


It  is  also  clear  that  (1  —  2ucpvp)2r  >  1  —  4ruicpi/p.  Combining  these  inequalities  gives 


(A. 11) 


(A. 12) 


bp  >  (1  -  4njcpi'p)(l  -  \CpVp) 


>  1  —  4(1  4-  ru)cpi/p.  I 


Proof  of  Lemma  A. 3 

Since  0  <  £p,r/p  <  1, 

~ £p  (1  ~  £p)  (1  —  ^p)  ~  Vp  (1  ~  £p)  (1  ~  Vp)  —  0- 

Factoring  the  lefthand  side  gives 

<fp7?p(2  —  <fp  —  Tip)  <  4p0  —  £p)  d"  7p(I  —  Tlp)- 


(A. 13) 


(A  14) 


(A. 15) 


Recall  that 


and 


"p- i 


Thus  by  (A. 15) 


4(2-{P-  VP) 

hl~i 


4(2  —  £p-i  —  Vp-i) 

*2,-. 


4(CP(1  zJp)  +  Vp(l-VP)) 

hl 


Vp-l 


< 


1 

‘fphp 


The  second  inequality,  (A. 6b),  is  clear  since 


v?  _  (1  ~  ZP)  4-  (1  -  rip) 

Vp—\  £p(  1  —  £p  )  4~  qp(  1  T)p) 

and  0  <  £p,  rjp  <  1.  I 


(A. 16) 


(A. 17) 
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Proof  of  Lemma  A. 4 


7  <  Cpi/p  <  T  <  1/4 


(A. 18a) 


47(1-27)  <  Cp — ! — 1  <  4r(l-2r). 


(A.  18b) 


Note  that  this  is  just  the  calculus  problem:  Find  the  maximum  and  minimum  of  /(x,  y)  = 
x(  1  —  x)  4-  y(l  —  y)  in  ft  =  {(x,  y):  27  <  x  +  y  <  2r,  x  >  0,  y  >  0}  ,  and  the  solution 
is  straightforward. 

By  induction,  it  is  easy  to  see  that 


^r+r  II  ( 1  ~  ^ 


^or  — n 


(A. 19) 


By  (A.  18)  this  is  true  for  n  =  1,  i.e. 


0  r,  W  ^  ^  0  „  2/3,  ^  0 

40,  V1  40+1)  -  cp-\Vp-l  ^  401  —  1  ^ 1  4a)  -  4o^T 


Assume  (A. 19)  is  true  for  cp_n+iup_n+i ,  then 


(A. 20) 


J=0  x 


4a  — n+2  XI  V  40+1-; 
>=0  v 


>_1_  jfA-JL 

—  4a— n+1  1  IX  y  4Q+1-J 


40—71+2 


(A. 21) 


40— n+i  xx  y  40+1— j  1 
j= 0  x  ' 


Using  []  U  ~  xj)  >  1  ~  E  gives 

J=0  ;=0 


nV  (e 


I* 


The  upper  bound  for  cp-nvp-n  is  obvious.  I 


Proof  of  Theorem  3.4  part  a. 

w(l  -w)  l2  ^  n{k)  ^  2ru?L2 
a')  max(2,d(l  +  rwj)^1  -  ~  'I-''1 ' 


for  v„  <  jy. 


From  lemmas  A. la  and  A. 2a,  p  =  k,  k  —  1, . . . ,  2, 


D\k)  <4,v  £c,  +X>i1). 


(A. 23) 


On  the  coarsest  grid,  D]  <  4rtuci .  Hence  by  the  definition  of  the  cp’s  (2.15) 


<  4rui  ^  c„  j  < 


16ru> 


(A. 24) 


Using  Ci  =  h\/ 8  gives  the  upper  bound 

D\k *  <  -ruhy 

To  get  the  lower  bound,  use  an  induction  argument.  By  lemma  A. 16, 


Let  1  <  p  <  k  and  assume  that 


D(p)  > 


u;(l  —  u)h\ 
max (2,  d(l  4-  ru))' 


(A.25) 


(A. 26) 


By  rearranging  terms,  using  lemmas  A. lb  and  A. 2b  (which  can  be  used  since  <  d/h\ 
implies  cp vp  <  d/( 8  ■  4P)  <  1/4  by  lemma  A. 4a)  it  is  seen  that: 


D<^'>  >  ■  _  +  4w(l  -  u.)c*.  (l  - 

-  max(2,  d(l  +rw))  1  j  P+1  V  max (2,  d(  1  4-  rw) 


;a.2S) 


Lemma  A. 3b  guarantees  that  up  <  Uk  <  jj  and  therefore  the  last  term  in  (A. 28)  is 

hx 

positive  and  can  be  thrown  out.  This  proves  part  a.)  of  Theorem  3.4.  I 


% 


Proof  of  Theorem  3.4  part  b. 

b.)  -Ml T. +  <  z,(‘>  <  i  for  „<*>  >  jl 

3(1  +  rui)vk  '  uk  h\ 

Using  the  definition  of  ck ,  vk  >  d/h\  implies  ckvk  >  d/( 2  ■  4fc+]).  For  each 
<  1 .  To  see  this,  first  note  that 


A;  =  1  -(1  -2wci£/1)2r  <  1. 


Lemma  A. 3a  together  with  the  definition  of  bp  imply 


(A. 29) 


(1  -  2ucpup)2rup_1 


Combining  (A. 29)  and  (A. 30)  with  the  definition  of  ap,  gives  D\p^  <\/up 
For  the  lower  bound,  divide  the  argument  into  two  cases.  Define 

7  =  [log4  2(1  +  ru)] 

where  [x]  means  the  greatest  integer  in  x. 


(A. 30) 


(A.31) 


case  1 


case  2 


CkVk  > 


2  •  ri' 


(A. 32) 


a 

—  <  ckvk  <  1  for  some  integer  a,  7  <  a  <  k.  (A. 33) 


By  the  definition  of  7, 


-(1  +  rw)  <  47  <  2(1  +  rw). 


(A. 34) 


For  case  i,  Lemma  A. lb  gives 


4w(l  —  w)  d  du>{  1  —  u>) 

ak  >  — 1 - L  >  —1 - -L.  (A. 35) 

vk  2  •  47  ( 1  +  ru)uk 

For  case  2,  look  at  the  finest  grid  on  which  the  eigenvector  is  “rough  enough"  For 
p  =  k  —  a  +  7,  lemma  A. 4  shows  that 


d  (  d  \  d 

2-47+1  3  •  r7+1  /  ~  CpVp  ~  2  4  7 


A.  36 ) 


V  V.V  v,  '  *  * 


'  y  y  ■ 


/  ^  .  J* 


Therefore  on  f lp , 


D\ ^  >  ap  >  4u;(l  —  <jj)cp  > 


d{  1  —  <jj)ui 


8(1  4-  ru)i/p 

Now  this  information  needs  to  get  back  to  the  fine  grid,  Q.k .  On  for  p  P 
A. 4  says 


1 


<  cnun  < 


2-4  or  — fc+p+i  3  .  ^a  —  k+p+l  J  —  P  P  —  2  •  40>-*  +  P  ' 

Now  using  lemmas  A.  16,  A. 26  and  A. 36  and  rearranging  terms, 


>  4w(l  -  u)cp  +  (1  -  4(1  +  ruj)cpi/p ) 


uj(  1  —  U>) 

8(1  +  rui)vp 


d(  1  — 

8(1  +  rui)vp 


4-  4w(l  —  u)cp 


1 


d' 

8 


> 


d(l  —  lj)uj 
8(1  4-  ru>)vp ' 


(A. 37) 
p,  lemma 

(A. 38) 

(A. 39) 


Since  this  is  true  for  all  p  >  p,  take  p  =  k.  I 
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